## script file for coin flip survey 
## 
require(plyr)
require(lattice)
require(grid)
require(wordcloud)
require(tm)
require(RColorBrewer)


s01 <- list()
s01 <- within(s01, {
    ## read in data, pull out locations of columns
    dat <- read.csv('s01.coins.csv')
    cols <- colnames(dat)
    ## get column pos for various responses
    c.hum1 <- grep('Coins.Human.1', cols)
    c.hum2 <- grep('Coins.Human.2', cols)
    c.coins <- grep('Coins.Random', cols)
    c.hurr <- grep('Hurricane', cols)
    ## time to complete survey
    timings <- data.frame(
        minutes= as.numeric(
            as.POSIXct(dat$Completed) - 
            as.POSIXct(dat$Date.started) 
        )
    )
})

## convenience plot -- visualize timings
plot.timings <- function(.dat, .logx=F) {
    densityplot(.dat,  
        plot.points='jitter', 
        xlab='Minutes to complete survey',
        scales=list(x=list(log=.logx)), 
        pch=20, cex=2, alpha=0.75
    )
}
    

s01.plots <- list()
s01.plots <- within(s01.plots, {
    ## how long to complete survey
    ## reverse order for lapply and plotting
    plot1.3 <- plot.timings(s01$timings$minutes, .logx=T) 
    plot1.2 <- plot.timings(subset(s01$timings, minutes<100)$minutes) 
    plot1.1 <- plot.timings(s01$timings$minutes) 
})





tally.trans <- function(vec) {
## get number of transitions and length of longest run
## for a single vector
    flips <- 0
    thisrun <- 0
    maxrun <- 0
    for(ii in 2:length(vec)) {
        ## not a flip, extend run
        if (vec[ii-1]  == vec[ii]) {
            thisrun <- thisrun + 1
        ## flip, update maxrun
        } else {
            maxrun <- max(maxrun, thisrun)
            flips <- flips + 1
            thisrun <- 1
        }
    } 
    return( data.frame( flips, maxrun))
}


## number of sims of 50 tosses
.nsim <- 1e3

s01 <- within(s01, {
    ## pull out hurricane replies as string vector
    dat.hurr <- unlist(dat[, c.hurr])
    
    ## pull out human and coins 
    dat.list <- list(
        Simulation = matrix( rbinom(.nsim*50, 1, 0.5), ncol=50),
        Humans = rbind( as.matrix(dat[, c.hum1]), as.matrix(dat[, c.hum2])),
        Coins = as.matrix(dat[, c.coins])
    )

    ## for each data.set
    dat.trans <- ldply(dat.list,  function(x) {
        ## for each row in dataset, tally trans and maxrun
        adply(x, 1, tally.trans)
    })
    ## relevel id so sims plot first, class data on top
    dat.trans$.id <- relevel(factor(dat.trans$.id), 'Simulation')
})
            

## max runs versus transitions for humans, coins, sims
s01.plots$plot2 <- 
    xyplot( jitter(maxrun) ~ jitter(flips), groups= .id, 
        s01$dat.trans, 
        cex=c(0.5, 2,2), pch=c(1,20,20), alpha=0.5, 
        auto.key=list(x=0.9, y=0.9, corner=c(1,1)), 
        xlab=sprintf('Number of transitions \n %s simulations', .nsim), ylab='Max Run Length'
)


pdf(file='s01.plots.pdf', width=6, height=4)
    lapply(s01.plots, function(x) print(x))
    ## wordcloud -- include all words, nice random colors
    wordcloud(s01$dat.hurr, min.freq=1, random.color=T, 
        rot.per=0.2, scale=c(2.5,0.35), 
        colors=brewer.pal(5, 'Set1')
    )
    grid.text('What factors are important in predicting hurricanes?', x=0.5,y=0.05, gp=gpar(fontsize=16))
dev.off()
